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I. Preface 


This document serves as a Final Report, summarizing the work performed from 
September 1, 1985 through January 30, 1989 under Grant NAG 3-645 from the 
NASA Lewis Research Center. The Principal Investigator for the Grant was Profes- 
sor David A. Caughey of the Sibley School of Mechanical and Aerospace Engineering 
of Cornell University; the Grant Technical Monitor was Dr. Rodrick V. Chima of 
the Computational Methods Branch of the NASA Lewis Research Center. 


II. Accomplishments 


Research was performed in the general area of Computational Aerodynamics, with 
particular emphasis on the development of efficient techniques for the solution of the 
Euler and Navier-Stokes equations for transonic flows through the complex blade 
passages associated with turbomachines. In particular, multigrid methods have been 
developed, using both explicit and implicit time-stepping schemes as smoothing 
algorithms. The specific accomplishments of the research supported by this grant 
have included: 
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(1) the development of an explicit multigrid method to solve the Euler equations 
for three-dimensional turbomachinery flows based upon the multigrid imple- 
mentation of Jameson’s explicit Runge-Kutta scheme [Jameson 1983]; 

(2) the development of an implicit multigrid scheme for the three-dimensional Euler 
equations based upon Lower-Upper factorization; 

(3) the development of a multigrid scheme using a diagonalized ADI Implicit algo- 
rithm; 

(4) the extension of the diagonalized ADI multigrid method to solve the Euler 
equations of inviscid flow for three-dimensional turbomachinery flows; and 

(5) the extension of the diagonalized ADI multigrid scheme to solve the Reynolds- 
averaged Navier-Stokes equations for two-dimensional turbomachinery flows. 

The work in each of these areas will now be described briefly. Copies of Smith & 
Caughey [1987], Yokota & Caughey [1988], Caughey [1988], and Caughey & Turkel 
I [1988] are also included as Appendices to this report, and contain further details 

regarding this work. 

The Grant has supported all or part of the Ph.D. Thesis research of Drs. Wayne 
A. Smith (Mechanical Enginering) and Jeffrey W. Yokota (Aerospace Engineering), 
who have completed their theses, and of Thomas Tysinger (Aerospace Engineering) 
and Culbert Laney (Applied Mathematics), both of whose thesis work is still in 
progress. 

I 


A. Explicit Runge-Kutta Multigrid Scheme 

A fast and efficient three dimensional Euler solver for transonic flow through rotating 
blade passages has been developed and tested. The time dependent equations are 
discretized spatially using a finite volume approach, and are advanced temporally 
with a multiple stage time stepping scheme, as popularized by Jameson, Schmidt, 
& Turkel [1981]; see also Jameson [1982]. A dramatic increase in the rate of con- 
vergence for steady state solutions has been achieved with a multigrid algorithm 
(Brandt [1973]) which employs the multistage scheme as its smoothing procedure. 
The ability of the multistage scheme to damp high frequency error is enhanced 
by a carefully selected combination of the coefficients governing implicit residual 
smoothing, added dissipation, time steps, and the multistage scheme itself. 

The Euler equations are formulated in terms of absolute flow variables in a 
Cartesian frame of reference rotating with constant angular velocity. This allows 
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the equations to be discretized in a way which ensures that a uniform stream satisfies 
the difference equations identically (Holmes & Tong [1985]). This property is not so 
important for the present turbomachinery application, but is critical for the solution 
of flows about propellers and helicopter blades. 

The use of centered differences ensures that the scheme is second order accurate 
provided the mesh is sufficiently smooth, but also allows the solution to decouple 
at odd- and even-numbered points in the grid, necessitating the addition of suitable 
dissipative terms. These are introduced as an adaptive blend of second- and fourth- 
differences which maintain high accuracy while capturing shock waves with little or 
no overshoot. 

A multiple stage time stepping scheme is used to discretize the remaining time 
derivative. The scheme is implemented within the context of the multigrid algorithm 
to accelerate convergence to the steady state. This is done following the procedure 
of Jameson [1983], but with the multigrid sequence defined recursively in a way 
which permits the definition of generalized “W-cycles,” and allows greater amounts 
of smoothing to be done effectively on the coarser grids of the sequence. 

The option to perform implicit residual smoothing was also included in the 
original code, but a numerical determination of the optimal values for use with 
multigrid of the stage coefficients and implicit smoothing parameter for a five-stage 
Runge-Kutta scheme showed that implicit smoothing should not be used (Smith & 
Caughey [1987], Smith [1987b]). 

A number of results have been computed, for both compressor and turbine ge- 
ometries; the values of the parameters determined from the simple one-dimensional 
optimization procedure have been used without modification to compute solutions 
to all the test cases. The grids used for these calculations were constructed using a 
modified version (Smith [1987a]) of the GRAPE code of Sorenson [1980]. 

Further details of the method and results of several calculations are included in 
Smith [1987b] and Smith & Caughey [1987]. These also include comparisons with 
the experimental and computational results of Chima & Strazisar [1982]. 


B. L-U Implicit Multigrid Scheme 

A Lower-Upper (L-U) Factored Implicit scheme for solving the Euler equations 
of inviscid flow in three-dimensional, rotating internal cascade geometries has also 
been developed. The Euler equations are written in conservation form and are 
transformed to generalized coordinates. These are approximated using a finite- 
volume discretization, and suitable artificial dissipative terms are added so that 
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shocks are captured automatically. The spatial discretization is similar to that of 
Jameson et al [1981] and to that described above. 

After a time-linearization of the nonlinear flux vectors, the implicit operator 
can be factored into lower- and upper- block bi-diagonal factors by using one-sided 
differences. This results in a scheme which requires the inversion only of the 5x5 
blocks at each cell for each of the two effectively explicit sweeps through the domain. 
This requires significantly less computational labor than the block tridiagonal inver- 
sions of conventional Alternating Direction Implicit schemes, and incorporates the 
added efficiency that only two sweeps are required even for three-dimensional prob- 
lems. The Jacobian matrices are split in a way which produces diagonally-dominant 
factors for each matrix inversion in each cell. This is ensured by constructing the 
splitting so that the eigenvalues of the forward-differenced matrices are non-positive, 
and those of the backward-differenced matrices are non-negative. The actual pro- 
cedure follows closely that suggested by Jameson & Turkel [1981]. 

Artificial dissipation terms are added to the centered scheme to reduce odd- 
even decoupling of the solution, and to suppress overshoots in the vicinity of shock 
waves; these terms are constructed following the formulation proposed by Jameson 
et al [1981]. The explicit terms are constructed of blends of nonlinear second- and 
fourth-differences, while only the second-differences are treated implicitly. 

The scheme is implemented within the context of the multigrid algorithm to 
accelerate convergence to the steady state. This is done following the procedure of 
Jameson [1983]. The calculations were performed on three-dimensional “H-type” 
grids for turbomachinery blade rows, which were generated using the version of the 
GRAPE code of Sorenson [1980] that was modified by Smith [1987a]. 

Further details of the method, as well as results illustrating its efficiency, are 
included in Yokota [1987] and Yokota & Caughey [1987]. 


C. Diagonal Alternating Direction Implicit Multigrid Scheme 

The development of Alternating Direction Implicit multigrid methods is an attrac- 
tive means to circumvent the problem of the slow convergence of explicit methods on 
grids containing cells of extremely high aspect ratio. These methods have also been 
studied recently by Jameson & Yoon [1985]. In order for the implicit method to be 
an effective smoothing algorithm when used in conjunction with the multigrid algo- 
rithm, it is important to include an accurate representation of the dissipative terms. 
Pulliam [1984] has also shown the importance of including the fourth-difference 
terms in the implicit operator even when multigrid is not used. Since, as described 
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in the preceding sections, these usually include fourth-differences to maintain high 
accuracy, their inclusion in the implicit operator requires the solution of block pen- 
tadiagonal systems for each one-dimensional factor. To avoid the high cost of solving 
block pentadiagonal systems, the equations in the present scheme are first diago- 
nalized at each point using a local similarity transformation, following Chaussee & 
Pulliam [1981]. This has the effect of decoupling the equations, and requiring the 
solution of four (in two-dimensional problems) scalar pentadiagonal equations for 
each factor. The resulting method has good high- wavenumber damping, so is a good 
smoothing algorithm for use in conjunction with the multigrid method. It is also 
computationally efficient because of the need to solve only scalar systems. Addi- 
tional computation is required to calculate the local similarity transformations, and 
to perform matrix multiplies of the residual, and of the intermediate and final cor- 
rections, but this is a small fraction of the work required to solve the block systems. 
In addition, on many computers these additional computations can be vectorized to 
further reduce the required CPU time. 

In preparation for the incorporation of the viscous terms in the multigrid ADI 
scheme for turbomachinery problems (see subsection II. E below), a version of the 
implicit multigrid Euler code for two-dimensional cascades was developed. This code 
has also been used to perform experiments demonstrating the importance of treating 
the dissipation in a fully conservative manner, for flows containing shock waves, even 
for internal flow problems. Results a,re summarized here for a test case involving 
moderately strong shock waves in the flow through a simple cascade of NACA 0012 
airfoils. The flow is computed through an unstaggered cascade of airfoils having 
a height-to-chord ratio of three for a ratio of downstream static to upstream total 
pressure of 0.631298 (corresponding tc isentropic expansion to a Mach number of 
Moo = 0.838). The grid for this calculation, generated using Smith’s extension 
(Smith [1987]) of the GRAPE Code (Sorenson [1980]), is shown in Figure 1. Blade 
surface distributions of surface pressure coefficient are shown in Figures 2(a) and 
2(b) for cases in which the dissipative terms were differenced non-conservatively and 
conservatively, respectively. Contours of constant pressure for these two cases are 
shown in Figures 3(a) and 3(b). As is evident from these figures, the shocks are 
significantly weakened and shifted upstream by the non-conservative differencing. 

The details of the new algorithm and results of several two-dimensional, tran- 
sonic airfoil calculations are described in Caughey [1987b, 1988]. In addition, new 
results concerning the effects of numerical dissipation in finite-volume schemes have 
been presented by Caughey [1987a] and by Caughey & Turkel [1988]. 
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D. Implicit Multigrid Scheme for the 3-D Euler Equations 

Work has also begun on development of an implicit multigrid method to solve the 
three dimensional Euler equations for the transonic flow through rotating blade 
passages. The time dependent equations are discretized spatially using a finite 
volume approximation (see Jameson et al [1981]), and are advanced temporally 
using an implicit, approximate-factorization time stepping scheme, as developed by 
Caughey [1988]. 

The Euler equations are formulated in terms of absolute flow variables in a 
Cartesian frame of reference rotating with constant angular velocity. This allows 
the equations to be discretized in a way which ensures that a uniform stream satisfies 
the difference equations identically (Holmes & Tong [1985]). This property is not so 
important for the present turbomachinery application, but is critical for the solution 
of flows about propellers and helicopter blades. 

The use of centered differences ensures that the scheme is second order accurate 
provided the mesh is sufficiently smooth, but also allows the solution to decouple 
at odd- and even-numbered points in the grid, necessitating the addition of suitable 
dissipative terms. These are introduced as an adaptive blend of second- and fourth- 
differences which maintain high accuracy while capturing shock waves with little or 
no overshoot. 

An implicit Alternating Direction Implicit (ADI) time stepping scheme is used 
to discretize the remaining time derivatives. Each of the factors of the ADI scheme 
is diagonalized by a local similarity transformation to decouple the linear equations 
to be solved for the corrections at each time step. This results in a scheme which is 
significantly more efficient than the usual form, which requires block pentadiagonal 
solutions for each of the three factors for each time step. The scheme is implemented 
within the context of the multigrid algorithm to accelerate convergence to the steady 
state. This is done following the procedure of Jameson [1983], but with the multigrid 
sequence defined recursively in a way which permits larger amounts of smoothing 
to be done on the coarser grids of the sequence. Under separate support, three- 
dimensional solutions for the steady aerodynamic problem of transonic flow past a 
swept wing have been developed using this method by Yadlin & Caughey [1988]. 

The implicit scheme has been implemented using the steady-state portions of 
the explicit multigrid code developed by Smith [1987b]. Further details of the steady 
state approximation and the multigrid method are described by Smith [1987b] and 
by Smith & Caughey [1987]. 

Preliminary results have been computed for a simple test geometry consisting 
of a symmetric channel between two NACA 0012 airfoil sections. The sections are 
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located three chord lengths apart, with no stagger. Plane walls representing the 
hub and shroud are located one chord apart. The grid used for these calculations 
was constructed using a version of the GRAPE code of Sorenson [1980], modified to 
allow construction of H-type grids (Smith [1987a]); separate two-dimensional grids 
were generated in each spanwise plane, then stacked to form the three-dimensional 
grid. Contours of constant pressure are shown in Figure 4 on the mid-channel 
blade-to-blade surface for a sample calculation. The downstream static to upstream 
total pressure was chosen to correspond to an inlet Mach number of approximately 
0.75. The calculation was performed on a grid containing 16 x 4 x 8 grid cells 
in the throughflow, radial, and blade-to-blade passage directions, respectively. A 
plot of the grid lines in the mid-channel plane is shown in the insert on Figure 4. 
Grids of higher density are, of course, necessary for meaningful results; this result 
is intended merely to show that the Diagonalized ADI multigrid solver is working 
for simple cases. The convergence history for this calculation is shown in Figure 
5. Plotted is the root - mean-square (over the entire domain) of the residual of the 
continuity equation, which reaches machine zero in approximately 160 Work units 
using two levels of multigrid. 


E. Implicit Multigrid Scheme for the 2-D Navier-Stokes Equations 

The multigrid diagonal implicit algorithm developed by Caughey [1988] for the Euler 
equations of inviscid, compressible flow, has been extended to solve the compress- 
ible Navier-Stokes equations in two-dimensions. Because of its stability properties, 
the implicit scheme is particularly attractive for computing solutions on grids with 
highly-stretched cells, as is required for high Reynolds number flows. Even for flows 
at moderate Reynolds number, it is sometimes necessary to use grids with high 
aspect ratio cells to resolve boundary layer regions. 

As for the Euler equation algorithm (Caughey [1988]), spatial derivatives are 
approximated using a finite volume formulation, and local time stepping is used to 
increase the convergence rate for steady problems. Artificial dissipation consisting of 
an adaptive blend of second and fourth differences is added to the scheme to insure 
convergence to a steady state and to allow accurate shock capturing for transonic 
flows. A recursive multigrid algorithm similar to that described by Smith & Caughey 
[1987] is implemented to accelerate convergence. 

In the implicit Euler algorithm, the time linearization of the contribution of 
the convective flux vectors to the implicit operator gives rise to Jacobians. These 
Jacobians can be diagonalized to increase the efficiency of the algorithm; this can 
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always be done by a similarity transformation, since the system describing inviscid 
flows is hyperbolic. For the Navier-Stokes equations, however, it is not possible both 
to include the viscous terms in the implicit factor and to diagonalize the system, 
since the convective and viscous Jacobians are not simultaneously diagonalizable. 
It is desirable to maintain the efficiency of a diagonalized scheme, and so an al- 
ternate solution is required. There has been some success by implicitly including 
an approximation to the viscous Jacobian eigenvalues (Pulliam [1986]), but for the 
present scheme the viscous terms are neglected completely in the implicit factor, 
and thus contribute only to the explicit part of the equation. This can have no 
effect on a converged solution; it can affect only the stability of the iteration and 
rates of convergence. Converged solutions have been obtained using such a scheme, 
although recent experience suggests that implicit treatment of the viscous terms is 
necessary to remove low-frequency oscillations in the solution which can develop in 
some cases. With diagonalization, the computational efficiency of the algorithm is 
comparable to that of the explicit multi-stage schemes - i.e., the CPU time per time 
step is virtually the same. 

The algorithm has been implemented in a computer code to calculate transonic 
flows past two-dimensional airfoils, and a number of test cases have been computed. 
The first case presented here is for the subsonic laminar flow (Re = 5000, M oo = 
0.5) past a two-dimensional NACA 0012 symmetric airfoil at zero degrees angle of 
attack. The calculation is performed on a 192 x 48 cell “C”-grid generated using 
the GRAPE code elliptic mesh generator (Sorenson [1980]). The outer boundary of 
the mesh is located about 3 chords from the body. Care is taken to insure sufficient 
clustering in the region close to the body surface where viscous effects are significant. 
Approximately 10 mesh points are included within the boundary layer at the airfoil 
trailing edge, and the first point normal tc the body surface is located at about 
.001 chords. The surface pressure distribution, presented in Figure 6, agrees well 
with that presented by Martinelli, Jameson, & Grasso [1986]. The flow separates 
at approximately 85% of the chord, as can be seen from the contour plot of the 
stream wise component of mass-flux density in Figure 7; this value is close to the 
values reported by both Swanson & Turkel [1985] and Jayaram & Jameson [1988] 
for this case. 

The iterative process is begun by initializing the solution to free stream values. 
A plot of the convergence history is shown in Figure 8. Using six levels of multigrid 
and local time stepping, the solution has converged to a steady state in approxi- 
mately 75 work units; this corresponds to 45 multigrid cycles. One work unit is 
defined as the amount of computational work required to advance the solution one 
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time step on the finest mesh level. It includes the work done on coarser levels during 
a multigrid cycle. Overall, the average residual is reduced by 7 orders of magnitude 
in about 300 work units or 180 multigrid cycles. This represents a substantial im- 
provement over rates obtained by other researchers (e.g., Martinelli et al [1986] and 
Swanson & Turkel [1985]. 

The second case described here is for the transonic laminar flow (Re = 
500, Moo = 0.8) past a NACA 0012 airfoil at ten degrees angle of incidence. A 
192 x 48 cell “C”-grid is again used, with its outer boundary located about 3 chords 
from the body. A spacing of .01 chords is used for the first grid line normal to the 
body surface. The surface pressure distribution and contours of constant streamwise 
mass-flux density are shown in Figures 9 and 10, respectively. The solution appears 
to be in reasonable agreement with that obtained by Martinelli, et al [1986]. In this 
case, the flow separates from the upper surface of the airfoil at approximately 35% 
of the chord, as can be seen in Figure 10, and a large recirculating region forms on 
the upper surface of the airfoil. 

The solution converges to a steady state from an initialized flow field of uniform 
free stream in less than 100 work units. The average residual is reduced by almost 10 
orders of magnitude in 300 work units as shown in Figure 11; this can be compared 
to the results reported by Martinelli et al [1986], where a convergence of only 5 
orders of magnitude in 1200 cycles v/as reported. 
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NACA 0012 Airfoil (C-Grid Cascade) 
Grid 128 x 32 


Figure 1. Computational grid for cascade test of non-conservative numerical dissi- 
pation. 


immimmuimmiiii iiiiiiiimmmmiiii 




Multigrid Calculation of Three-Dimensional Turbomachinery Flows 15 



Figure 2(a). Blade surface pressure distribution for transonic flow through cascade; 
non-conservative numerical dissipation. 



Figure 2(b). Blade surface pressure distribution for transonic flow through cascade; 
fully conservative numerical dissipation. 
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Figure 4. Contours of constant pressure at mid-plane of three-dimensional test 
channel. 
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Figure 5. Residual convergence history for three-dimensional diagonal ADI multigrid 
solution. 
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FL053C: NACA 0012 FREEAIR SUBSONIC 
Mach 0.500 Alpha 0.000 Re 5000.0 

Cl 0.0000 Cd 0.0495 Cm 0.0000 

Grid 192x48 Work 300.88 Res0.601E-10 


Figure 6. Surface pressure distribution for symmetrical, laminar flow past NACA 
0012 Airfoil; Re = 5000.0, = 0.50, a = 0.0 
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FL053C: NACA 0012 FREEAIR SUBSONIC 


Minimum = 1000E+00 Incrmnt - 0.2500E-01 

Maximum = 0.7000E+00 W(~,i>j) COntOUTS Scale =0.1 OOOE+05 

Figure 7. Streamwise mass-flux density for symmetrical, laminar flow past NACA 
0012 Airfoil; Re = 5000.0, Moo = 0.50, a = 0.0 
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Figure 8. Residual convergence history for symmetrical, laminar flow past NACA 
0012 Airfoil; Re = 5000.0, = 0.50, cc = 0.0 
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FL053C: NACA 0012 FREEAIR T.L.F. 

Mach 0.800 Alpha 10.000 Re 500.0 

Cl 0.4277 Cd 0.2407 Cm 0.0000 

Grid 192x48 Work 300.88 Res0.517E-12 


Figure 9. Surface pressure distribution for transonic, separated flow past NACA 
0012 Airfoil; Re = 500.0, Moo = 0.80, a = 10.0 
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FL053C: NACA 0012 FREEAIR T.L.F. 


Minimum = -.2000E+00 Incrmnt = 0.1 000E+00 

Maximum = 0.1000E+01 w(2,i,j) COIltOUrS Scale =0.1000E+01 

Figure 10. Streamwise mass-flux density for transonic, separated flow past NACA 
0012 Airfoil; Re — 500.0, M«x, = 0.80, a = 10.0 
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Mach 0.800 Alpha 10.000 

Resl 0.120E-02 Re 500. 
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CFL 4.00 
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Figure 11. Residual convergence history for transonic, separated flow past NACA 
0012 Airfoil; Re = 500.0, M ^ = 0.80, a = 10.0 




Multigrid Solution of the Euler Equations 
Wayne Allen Smith, Ph.D. 

Cornell University 1987 


A fast solver for three-dimensional, inviscid, transonic flow in rotating 
domains is presented. The scheme is designed to be to be extended to the 
numerical simulation of viscous transonic flows about realistic geometries. 
The techniques employed are therefore chosen for their extensibility as well 
as for their efficiency. Techniques that rely on simplifications made to more 
general flow equations by assumptions of inviscid flow of an ideal gas are ruled 
out, as are those whose analytic manipulations of the governing equations 
would become unwieldy as the complexity of the equations increases. Even 
with this generality, the resulting scheme provides the most efficient solution 
of the Euler equations for internal flows available to date. 

Spatial discretization of the governing equations is performed with a finite 
volume scheme to provide geometric generality and to facilitate the construc- 
tion of a conservative scheme. Temporal discretization is performed with a 
multiple-stage time-stepping scheme to allow the calculation of both time- 
accurate and steady solutions. The use of multigrid plays a critical role 
in making the explicit scheme efficient. The effects of residual smoothing 
for convergence acceleration are examined and residual smoothing is found 
to be detrimental when used in conjunction with multigrid and an efficient 
high wavenumber damping multistage scheme. Multistage parameters that 


provide efficient high wavenumber schemes are determined through an opti- 
mization procedure that uses a simple hyperbolic model equation. Boundary 
conditions employing simple extrapolations based solely on assumptions of 
the hyperbolic character of the governing equations are shown to be suffi- 
cient. The efficiency of the scheme is demonstrated with representative com- 
putational results taken from a high speed compressor. Finally, the great 
potential of the scheme to exploit computational parallelism is demonstrated 
by implementing the scheme on a hypercube, a local-memory multi-processor 
computer. 


